#Alex Gazmararian
#afg2@princeton.edu
#January 18, 2024

#Load packages
library(tidyverse)
library(modelsummary)
library(kableExtra)
library(haven)
library(survey)
library(here)
library(viridis)

#Load function for custom latex
source(here("code", "fun", "fix_txt.R"))

#Load data
g <- readRDS(here("data", "CCES", "ces_processed.rds"))
#Create survey object, which will be used later for estimating weighted means for some variables
g.svy <- svydesign(ids=~caseid,weights=~teamweight,data=g)
#Setup table gof statistics
gm <- gof_map
gm[gm$raw=="nobs",]$clean <- "$N$"
gm[gm$raw=="adj.r.squared",]$clean <- "Adjusted $R^2$"
gm[gm$clean=="R2",]$omit <- TRUE
gm[c(16:70),]$omit <- TRUE

#Create Table 1: willingness to pay for solar
g %>%
  filter(!is.na(wtp)) %>%
  group_by(wtp,solarinterest) %>%
  summarize(n=sum(weight)) %>%
  group_by(solarinterest) %>%
  mutate(pct = prop.table(n) * 100,
         solarinterest = ifelse(solarinterest == 1, "Yes", "No")) %>%
  dplyr::select(-n) %>%
  pivot_wider(names_from = solarinterest, values_from = pct) %>%
  kbl(digits = 0, col.names = c(" ", "No", "Yes"), booktabs = TRUE, linesep = "", format = "latex", escape = FALSE) %>%
  add_header_above(c(" " = 1, "Solar Interest" = 2)) %>%
  cat(., file = here("output", "tables", "tab_pay_agree.txt"))
#Chi-square test for table note
chisq.test(g$wtp,g$solarinterest)

#Check how many would donate among solar interested population
g %>% filter(solarinterest==1) %>%
  mutate(yesdonate = ifelse(self_donate>0,1,0)) %>%
  group_by(yesdonate) %>%
  summarize(n = sum(teamweight)) %>%
  mutate(prop=n/sum(n))

#Figure 3: Allocation of net metering proceeds
plot.donate <- g %>%
  filter(!is.na(self_donate)) %>%
  pivot_longer(cols=c(self_donate,others_donate)) %>%
  mutate(
    name=ifelse(name=="others_donate", "Neighbors", "Own"),
    name=factor(name,ordered = TRUE, levels = c("Own","Neighbors")),
    value=(value/20)*100,
    solarinterest = ifelse(solarinterest==1,"Solar Interest","No Solar Interest")
  ) %>%
  ggplot(aes(x=name,y=value,fill=solarinterest,weight=weight)) +
  geom_violin() +
  facet_wrap(~solarinterest) +
  stat_summary(fun.y=mean, geom="point", size=3, aes(group=solarinterest,shape=solarinterest),
               fill="white",
               position = position_dodge(.9)) +
  scale_shape_manual(values=c(21,21)) +
  stat_summary(aes(label=paste0(round(..y..,0),"%")),
               fun.y=mean, geom="text", size=5, position=position_dodge(.9),
               hjust = -1) +
  theme_bw(base_size = 14) +
  viridis::scale_fill_viridis(discrete=TRUE) +
  scale_y_continuous(labels = scales::percent_format(scale=1), expand = c(0,0)) +
  labs(x="Expectations of Own versus Neighbors' Donation Allocation",
       y="Share of $20 Donated (weighted)",fill="Interest\nin Solar") +
  guides(shape=FALSE) +
  theme(
    legend.position = "",
    panel.grid = element_blank()
  )
plot.donate
ggsave(
  plot.donate,
  filename = here("output", "figures", "donate.png"),
  scale = 1.5,
  dpi = 300,
  width = 6.5,
  height = 3
)

#Create Table 2: Regression of donation allocation on individual covariates
ols.donate <- list()
ols.donate[["(1)"]] <- lm(
    scale(self_donate) ~ woman + scale(age) + income5 + race4 + partysummary + college + rural,
    weights = weight,
    data = g
    )
ols.donate[["(2)"]] <- lm(
  scale(self_donate) ~ woman + scale(age) + income5 + race4 + partysummary + college + gw_index + rural,
  weights = weight,
  data = g
)
ols.donate[["(3)"]] <- lm(
  scale(self_donate) ~ woman + scale(age) + income5 + race4 + partysummary + college + solarinterest + rural,
  weights = weight,
  data = g
)
ols.donate[["(4)"]] <- lm(
  scale(self_donate) ~ woman + scale(age) + income5 + race4 + partysummary + college + solar_inequity + rural,
  weights = weight,
  data = g
)
ols.donate[["(5)"]] <- lm(
  scale(self_donate) ~ woman + scale(age) + income5 + race4 + partysummary + college + gw_index + solarinterest + solar_inequity + rural,
  weights = weight,
  data = g
)
#estimate weighted donation average and standard deviation for table bottom
mean.outcome <- round(survey::svymean(~self_donate, g.svy, na.rm=TRUE)[[1]],3)
std.outcome <- round(sqrt(survey::svyvar(~self_donate, g.svy, na.rm=TRUE))[[1]],3)
#Create table
modelsummary(
  ols.donate,
  stars = c("*"=.1,"**"=.05,"***"=.01),
  vcov = "HC3",
  coef_map = c(
    "gw_index"="Climate change concern (index)",
    "solarinterest"="Intend to buy solar",
    "solar_inequity"="Rich benefit from solar",
    "woman"="Woman",
    "scale(age)"="Age (standardized)",
    "college"="College Degree",
    "rural"="Rural",
    "race4Black"="Black",
    "race4Hispanic"="Hispanic",
    "race4Other"="Other",
    "income5Less than $29,999"="Less than \\$29,999",
    "income5$30,000 - $59,999"="\\$30,000-59,999",
    "income5$100,000 or more"="\\$100,000 or more",
    "income5Prefer not to say"="Prefer not to say",
    "partysummaryNeither"="Neither",
    "partysummaryRepublican"="Republican",
    "(Intercept)"="(Intercept)"
  ),
  gof_map = gm,
  add_rows = data.frame(
    c("Donation Average", "Donation Standard Deviation"),
    c(mean.outcome, std.outcome),
    c(mean.outcome, std.outcome),
    c(mean.outcome, std.outcome),
    c(mean.outcome, std.outcome),
    c(mean.outcome, std.outcome)
  ),
  #output = "latex",
  escape = FALSE,
  title = "Regression of Donation to Help Install and Maintain Community Solar in Marginalized Communities on Individual Covariates in a Nationally Representative Sample of the American Public \\label{tab:donate}",
  notes = list("Notes: HC3 standard errors reported. Regressions employ survey weights.")
) %>%
  kableExtra::add_header_above(c(" "= 1,"DV: Donation to Install Solar (standardized)" = 5)) %>%
  kableExtra::pack_rows("Race (Baseline: White)", 15, 20) %>%
  kableExtra::pack_rows("Income (Baseline: \\$60,000-99,999)", 21, 28, escape = FALSE) %>%
  kableExtra::pack_rows("Party (Baseline: Democrat)", 29, 31) %>%
  kableExtra::kable_styling(font = 10) %>%
  cat(., file = here("output", "tables", "ols_donation.txt"))

#Table B1: Robustness of results when using ideology as a covariate  

ols.ideo <- list()
ols.ideo[["(1)"]] <- lm(
  scale(self_donate) ~ woman + scale(age) + income5 + race4 + ideo3 + college + rural,
  weights = weight,
  data = g
)
ols.ideo[["(2)"]] <- lm(
  scale(self_donate) ~ woman + scale(age) + income5 + race4 + ideo3 + college + gw_index + rural,
  weights = weight,
  data = g
)
ols.ideo[["(3)"]] <- lm(
  scale(self_donate) ~ woman + scale(age) + income5 + race4 + ideo3 + college + solarinterest + rural,
  weights = weight,
  data = g
)
ols.ideo[["(4)"]] <- lm(
  scale(self_donate) ~ woman + scale(age) + income5 + race4 + ideo3 + college + solar_inequity + rural,
  weights = weight,
  data = g
)
ols.ideo[["(5)"]] <- lm(
  scale(self_donate) ~ woman + scale(age) + income5 + race4 + ideo3 + college + gw_index + solarinterest + solar_inequity + rural,
  weights = weight,
  data = g
)
#Specify file path to save
tab_ideo <- here("output", "tables", "ols_donation.ideo.txt")
#Create table
modelsummary(
  ols.ideo,
  stars = c("*"=.1,"**"=.05,"***"=.01),
  vcov = "HC3",
  coef_map = c(
    "(Intercept)"="Intercept",
    "gw_index"="Climate change concern (index)",
    "solarinterest"="Intend to buy solar",
    "solar_inequity"="Rich benefit from solar",
    "woman"="Woman",
    "scale(age)"="Age (standardized)",
    "college"="College Degree",
    "rural"="Rural",
    "race4Black"="Black",
    "race4Hispanic"="Hispanic",
    "race4Other"="Other",
    "income5Less than $29,999"="Less than \\$29,999",
    "income5$30,000 - $59,999"="\\$30,000-59,999",
    "income5$100,000 or more"="\\$100,000 or more",
    "income5Prefer not to say"="Prefer not to say",
    "ideo3Neither"="Neither ideology",
    "ideo3Conservative"="Conservative"
  ),
  gof_map = gm,
  add_rows = data.frame(
    c("Donation Average", "Donation Standard Deviation"),
    c(mean.outcome, std.outcome),
    c(mean.outcome, std.outcome),
    c(mean.outcome, std.outcome),
    c(mean.outcome, std.outcome),
    c(mean.outcome, std.outcome)
  ),
  fmt = 2,
  output = "latex",
  escape = FALSE,
) %>%
  add_header_above(c(" "= 1,"DV: Donation to Install Solar (standardized)" = 5)) %>%
  cat(., file = tab_ideo)
#Format latex table
fix_txt(tab_ideo)

#Table B2: Robustness to using a control for both ideology and partisan identification
ols.ideopid <- list()
ols.ideopid[["(1)"]] <- lm(
  scale(self_donate) ~ woman + scale(age) + income5 + race4 + ideo3 + partysummary + college + rural,
  weights = weight,
  data = g
)
ols.ideopid[["(2)"]] <- lm(
  scale(self_donate) ~ woman + scale(age) + income5 + race4 + ideo3 + partysummary + college + gw_index + rural,
  weights = weight,
  data = g
)
ols.ideopid[["(3)"]] <- lm(
  scale(self_donate) ~ woman + scale(age) + income5 + race4 + ideo3 + partysummary + college + solarinterest + rural,
  weights = weight,
  data = g
)
ols.ideopid[["(4)"]] <- lm(
  scale(self_donate) ~ woman + scale(age) + income5 + race4 + ideo3 + partysummary + college + solar_inequity + rural,
  weights = weight,
  data = g
)
ols.ideopid[["(5)"]] <- lm(
  scale(self_donate) ~ woman + scale(age) + income5 + race4 + ideo3 + partysummary + college + gw_index + solarinterest + solar_inequity + rural,
  weights = weight,
  data = g
)
#create table
tab_ideopid <- here("output", "tables", "ols_donation_ideopluspid.txt")
modelsummary(
  ols.ideopid,
  stars = c("*"=.1,"**"=.05,"***"=.01),
  vcov = "HC3",
  coef_map = c(
    "(Intercept)"="Intercept",
    "gw_index"="Climate change concern (index)",
    "solarinterest"="Intend to buy solar",
    "solar_inequity"="Rich benefit from solar",
    "woman"="Woman",
    "scale(age)"="Age (standardized)",
    "college"="College Degree",
    "rural"="Rural",
    "race4Black"="Black",
    "race4Hispanic"="Hispanic",
    "race4Other"="Other",
    "income5Less than $29,999"="Less than \\$29,999",
    "income5$30,000 - $59,999"="\\$30,000-59,999",
    "income5$100,000 or more"="\\$100,000 or more",
    "income5Prefer not to say"="Prefer not to say",
    "ideo3Neither"="Neither Ideology",
    "ideo3Conservative"="Conservative",
    "partysummaryNeither"="Neither Party",
    "partysummaryRepublican"="Republican"
  ),
  gof_map = gm,
  add_rows = data.frame(
    c("Donation Average", "Donation Standard Deviation"),
    c(mean.outcome, std.outcome),
    c(mean.outcome, std.outcome),
    c(mean.outcome, std.outcome),
    c(mean.outcome, std.outcome),
    c(mean.outcome, std.outcome)
  ),
  fmt = 2,
  output = "latex",
  escape = FALSE,
) %>%
  add_header_above(c(" "= 1,"DV: Donation to Install Solar (standardized)" = 5)) %>%
  cat(., file = tab_ideopid)
fix_txt(tab_ideopid)

#Table B3: Robustness for a joint measure of ideology and partisanship
ols.joint <- list()
ols.joint[["(1)"]] <- lm(
  scale(self_donate) ~ woman + scale(age) + income5 + race4 + ideo_pid + college + rural,
  weights = weight,
  data = g
)
ols.joint[["(2)"]] <- lm(
  scale(self_donate) ~ woman + scale(age) + income5 + race4 + ideo_pid + college + gw_index + rural,
  weights = weight,
  data = g
)
ols.joint[["(3)"]] <- lm(
  scale(self_donate) ~ woman + scale(age) + income5 + race4 + ideo_pid + college + solarinterest + rural,
  weights = weight,
  data = g
)
ols.joint[["(4)"]] <- lm(
  scale(self_donate) ~ woman + scale(age) + income5 + race4 + ideo_pid + college + solar_inequity + rural,
  weights = weight,
  data = g
)
ols.joint[["(5)"]] <- lm(
  scale(self_donate) ~ woman + scale(age) + income5 + race4 + ideo_pid + college + gw_index + solarinterest + solar_inequity + rural,
  weights = weight,
  data = g
)
#create table
tab_joint <- here("output", "tables", "ols_donation_joint.txt")
modelsummary(
  ols.joint,
  stars = c("*"=.1,"**"=.05,"***"=.01),
  vcov = "HC3",
  coef_map = c(
    "(Intercept)"="Intercept",
    "gw_index"="Climate change concern (index)",
    "solarinterest"="Intend to buy solar",
    "solar_inequity"="Rich benefit from solar",
    "woman"="Woman",
    "scale(age)"="Age (standardized)",
    "college"="College Degree",
    "rural"="Rural",
    "race4Black"="Black",
    "race4Hispanic"="Hispanic",
    "race4Other"="Other",
    "income5Less than $29,999"="Less than \\$29,999",
    "income5$30,000 - $59,999"="\\$30,000-59,999",
    "income5$100,000 or more"="\\$100,000 or more",
    "income5Prefer not to say"="Prefer not to say",
    "ideo_pidConservative Democrat" = "Conservative Democrat",
    "ideo_pidConservative Independent" = "Conservative Independent",
    "ideo_pidConservative Republican" = "Conservative Republican",
    "ideo_pidLiberal Independent" = "Liberal Independent",
    "ideo_pidLiberal Republican" = "Liberal Republican",
    "ideo_pidModerate Democrat" = "Moderate Democrat",
    "ideo_pidModerate Independent" = "Moderate Independent",
    "ideo_pidModerate Republican" = "Moderate Republican"
  ),
  gof_map = gm,
  add_rows = data.frame(
    c("Donation Average", "Donation Standard Deviation"),
    c(mean.outcome, std.outcome),
    c(mean.outcome, std.outcome),
    c(mean.outcome, std.outcome),
    c(mean.outcome, std.outcome),
    c(mean.outcome, std.outcome)
  ),
  fmt = 2,
  output = "latex",
  escape = FALSE,
) %>%
  add_header_above(c(" "= 1,"DV: Donation to Install Solar (standardized)" = 5)) %>%
  cat(., file = tab_joint)
fix_txt(tab_joint)

#Figure C1: American intentions to electrify their lives
df.buy <- g %>%
  tidyr::pivot_longer(cols=c(home_stove:home_hybrid)) %>%
  dplyr::filter(!is.na(value)) %>%
  dplyr::group_by(name,value) %>%
  dplyr::summarise(n=sum(weight)) %>%
  dplyr::group_by(name) %>%
  dplyr::mutate(pct=n/sum(n),
                value = ifelse(value=="I don't know what this is", "Don't know", as.character(value)),
                value=factor(value,ordered=TRUE,levels=c("Already have","Yes","No","Don't know")),
                name=dplyr::case_when(
                  name=="home_waterheater"~"On-demand electric hot water",
                  name=="home_stove"~"Induction stove",
                  name=="home_solar"~"Solar panels",
                  name=="home_hybrid"~"Hybrid Car",
                  name=="home_heatpump"~"Heat pump",
                  name=="home_ev"~"All Electric Car"
                ))
df.buy.order <- df.buy %>%
  dplyr::filter(value=="Yes") %>%
  dplyr::arrange(desc(pct))
plot.buy <-
  df.buy %>%
  dplyr::mutate(pct_label = pct,
                pct_label = ifelse(pct_label<0.02, NA, pct_label)) %>%
  dplyr::mutate(name=factor(name,ordered=TRUE,levels=rev(df.buy.order$name))) %>%
  ggplot(aes(x=name,y=pct,fill=value,label=scales::percent(pct_label,accuracy=1))) +
  geom_col() +
  coord_flip() +
  theme_classic(base_size=14) +
  scale_y_continuous(labels=scales::percent, expand=c(0,0)) +
  viridis::scale_fill_viridis(discrete=TRUE,option="D") +
  geom_text(size = 5, fontface = "bold", aes(color=value),
            position = position_stack(vjust = 0.5)) +
  scale_color_manual(values=c("white","white","black","black")) +
  guides(color=FALSE) +
  labs(color="",y="Percent (Weighted)",x="Product",fill=stringr::str_wrap("Considering",20))
plot.buy
ggsave(
  plot.buy,
  filename = here("output", "figures", "electrifyintentions.png"),
  scale = 1.5,
  dpi = 300,
  width = 6.5,
  height = 3
)

#Figure 2: beliefs about which societal groups have benefited from or been harmed by solar deployment
df.dist <- g %>%
  tidyr::pivot_longer(cols=c(solar_poor:solar_firms)) %>%
  dplyr::filter(!is.na(value)) %>%
  dplyr::mutate(
    name = dplyr::case_when(
      name=="solar_urban"~"Those living in urban areas",
      name=="solar_rural"~"Those living in rural areas",
      name=="solar_rich"~"High-income earners",
      name=="solar_poor"~"Low-income earners",
      name=="solar_poc"~"Racial and ethnic minorities",
      name=="solar_median"~"The average person",
      name=="solar_firms"~"Energy companies"
    ),
    value = dplyr::case_when(
      value%in%c("Harmed a lot", "Harmed a little")~"Harmed",
      value%in%c("Benefited a little","Benefited a lot")~"Benefited",
      T~as.character(value)
    )
  ) %>%
  dplyr::group_by(name,value) %>%
  dplyr::summarize(n=sum(weight)) %>%
  dplyr::group_by(name) %>%
  dplyr::mutate(pct=n/sum(n))
#create vector with order of who has benefited
dist.order <-
  df.dist %>%
  dplyr::filter(value=="Benefited") %>%
  dplyr::arrange(desc(pct))
plot.dist <-
  df.dist %>%
  dplyr::mutate(
    name = factor(name,ordered=TRUE,levels=rev(dist.order$name)),
    value = factor(value,ordered=TRUE,levels=rev(c("Harmed","No effect","Benefited")))
  ) %>%
  ggplot(aes(x=name,y=pct,fill=value,label=scales::percent(pct,accuracy=1))) +
  geom_col() +
  coord_flip() +
  theme_classic(base_size=14) +
  scale_y_continuous(labels=scales::percent,expand=c(0,0)) +
  viridis::scale_fill_viridis(discrete=TRUE,option="E") +
  geom_text(size = 4, fontface = "bold", aes(color=value),
            position = position_stack(vjust = 0.5)) +
  scale_color_manual(values=c("white", "black","black")) +
  labs(color="",y="Percent (Weighted)",x="Group",fill="") +
  guides(color=FALSE)
plot.dist
ggsave(
  plot.dist,
  filename = here("output", "figures", "solarbenefit.png"),
  scale = 1.5,
  dpi = 300,
  width = 6.5,
  height = 3
)
